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1. INTRODUCTION POOR QUALITY 

The utility of numerical methods for predicting transonic flows 
* over' v;ings and bodies is by now well established. The computer pro- 
gram FL022, based on a method presented at the 1973 IFIP Symposium on 
Computing Methods [1] / has actually been v/idely used to calculate the 
aerodynaraic performance of v/ings of transport aircraft. Provided that 
a correction is made for the displacement effect of the viscous bound- 
ary layer, this code has been found to give predictions which are 
accurate enough to serve as a useful design guide [ 2 ] . The salient 
features of the code are: 

(1) the use of a potential flow approximation to the equations of 
motion 

(2) the use of upwind differencing in the supersonic zone to simulate 
the region of dependence of the flow, and to prevent the appear- 
ance of expansion shock waves which would violate the entropy 
inequality 

♦ 

(3) the use of a relaxation procedure based on an artificial time 
dependent equation to solve the difference equations 

(4) the use of a curvilinear coordinate system generated by a sequence 
of simple transformations to -produce coordinate surfaces following 
the wing shape. 

The use of the potential flov; approximation greatly reduces the 
amount of computation required. Since the resulting flow is irrotation- 
al, it is consistent to approximate shock waves by discontinuities 
across v/hich entropy is conserved. This approximation has been found 
quite satisfactory in practice, since the shock vjaves generated by air- 
planes cruising at subsonic speeds are generally quite weak. In fact 
the appearance of stronger shock v/aves marks the onset of drag rise , 
which sets an upper bound on the cruising speed . In order to obtain a 
unique solution to the potential flow equation, it is necessary to 

* This work v/as supported by the Office of Naval Research under Contract 
N00014-77-C-0032, and also by NASA under Grants NGR 33-016-167 and 
NGR 33-016-201. The calculations v/ere performed at the ERDA Mathe- 
matics and Computing Laboratory, under Contract EY-76-C-02-3077 . *000 . 
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1. INTRODUCTION QUaUTY 

The utility of nurrerical methods for predicting transonic flows 
• over wings and bodies is by now well established. The computer pro- 
gram FL022, based on a method presented at the 1973 IFIP Symposium on 
Coniputing Methods (I), has actually been v;idely used to calculate the 
aerodynamic performance of wings of transport aircraft. Provided that 
a correction is raade for the displacement effect of the viscous bound- 
ary layer, this code has been found to give predictions which are 
accurate enough to serve as a useful design guide (2]. The salient 
features of the code are: 

(1) the use of a potential flow approximation to the equations of 
motion 

(2) the use of upwind differencing in the supersonic zone to simulate 
the region of dependence of the flow, and to prevent the appear- 
ance of expansion shock waves which would violate the entropy 
inequality 

(3) the use of a relaxation procedure based on an artificial time 
dependent equation to solve the difference equations 

(4) the use of a curvilinear coordinate system generated by a sequence 
of simple transformations to produce coordinate surfaces following 
the wing shape. 

The use of the potential flov; approximation greatly reduces the 
eunount of computation required. Since the resulting flov: is irrotation- 
al, it is consistent to approximate shock waves by discontinuities 
across v;hich entropy is conserved. This approximation has been found 
quite satisfactory in practice, since the shock waves generated by air- 
planes cruising at subsonic speeds are generally quite weak. In fact 
the appearance of stronger shock v;aves marks the onset of drag rise, 
which sots an upper bound on the cruising speed. In order to obtain a 
unique solution to the potential flow equation, it is necessary to 

* This work v/as supported by the Office of Naval Research under Contract 
N00014-77-C-0032 , and also by NASA under Grants NGR 33-016-167 and 
NGR 33-016-201 . The calculations v/ere performed at the 2RDA Mathe- 
matics and Computing Laboratory, under Contract EY-76-C-02-3077 . *000 . 


exclude expansion shock v/avcs, corresiponding to the condition that 
entropy can only increase. The use of upwind differencing in the 
supersonic zone, first introduced by Murman and Cole [3] , has been 
found an effective way to enforce the entropy condition. The non- 
linear equeitions generated by the discrete approximation are not easy 
to solve. The use of a relaxation process modeled on an artificial 
time dependent equation [4] has been found to give reliable and 
acceptably fast convergence. 

The main disadvantages of the scheme used in FL022 are the use of 
nonconservative difference formulas, which result in a failure, to sat- 
isfy conservation of mass across shock waves, and the difficulty of 
finding suitable transformations of coordinates to permit the treatment 
of more complex geometric configurations. The method to be described 
here is an attempt to overcome these shortcomings, while retaining the 
successful features of the previous method. The basic idea is to use 
a discrete approximation which directly represents a balance of the 
mass flow through small volume elements. This leads to a relatively 
simple treatment of the potential flow equation in conservation form. 

The volume elements arc distorted cubes generated by local trilinear 
transformations defined by the element vertices. Elements of this kind 
can be packed around any reasonably smooth configuration. The subsonic 
difference formulas can conveniently be derived from the Bateman varia- 
tional principle [5] . A directional bias is introduced in the super- 
sonic zone by adding an artificial viscosity, which is constructed in 
such a way as to produce an effective switch to upwind differencing. 

This serves to prevent the appearance of expansion shock waves. The 
artificial viscosity has a divergence form, so chat the conservation 
form of the equations is preserved by the difference scheme, and proper 
shock jump relations, consistent v/ith the isentropic approximation, 
are satisfied in the limit as the mesh v?idth is decreased to zero [ 6 ] . 
The most promising alternative to the use of artificial viscosity to 
enforce the entropy condition appears to be the optimal control method 
proposed by Glowinski and Pironneau [ 1 \ , in which the entropy condition 
is represented by penalty functions . 

CVi?R)vNAI. PAi.ii’- ill 

2. FORMULATION OF THE EQUATIONS OF POOR QUALUT 

The flow is assumed to be isentropic and to satisfy the equations 
of potential flow. Let q be the velocity vector, with magnitude q , 
and P the density. Then the potential flow equation can be written in 
conservation form as 


( 1 ) 


V* (pq) “ 0 


where q is the gradient o£ the p< Lcntial 


( 2 ) 


q = Vcf. 


Let a be the local. speed of sound, and M the Kach number q/a. Also lot 
, qc^ = 1 and = 1 be the Mach nxmber, speed and density of the 
uniform flow at infinity. Then the local density is given by the 
formula 


( v-1 2/- 2^ 1 

(3) p = |l + H‘(l - q^)| 

t 

where, y is the ratio of specific heats, and the pressure and speed of 
sound follow from the relations 


(4) 



Y-1 



Equation (1) is hyperbolic in supersonic flow {M > 1) and elliptic 
in subsonic flow, and shock waves v;ill generally appear if there is a 
region of supersonic flow. The shock jump conditions are 

(a) continuity of i|) implying continuity of the tangential velocity 
component 

(b) continuity of pq , where g is the normal velocity component 

(c) the entropy condition that g^ decreases through the shock. 

Under the assumption of isentropic flov;, conditions (a) and (b) imply 
that the normal component of momentum is not conserved. The resuiting 
momentum deficiency causes the appearance of a drag force, v/hich is 
an approximat-ion to the wave drag [8] . 

The boundary condition at the body is 


(5) 



To obtain a unique lifting solution we also impose the Kutta condition 
that the flov; leaves the trailing edge smoothly with equal velocities 
along the upper and lower surfaces. The resulting spanwise variation 
in the circulation T = q ds around each section of the wing causes 
a vortex sheet to be shed from the trailing edge. The vortex sheet 
will be convected v;ith the flov7,and roll up along its side edges. In 
the calculations this will be ignored and the vortex sheet will be 
assumed to coincide with a coordinate surfare. The conditions applied 


at the sheet are then 

(a) the jump V in the potential is constant along lines parallel to 
the free stream 

(b) the normal velocity component is continuous through the shock. 
According to an analysis of the asymptotic behavior of the potential 
in the far field [9] , <}) approaches the potential of the undisturbed 
uniforra flow except in the Trefftz plane far dovmstream, where it 
satisfies the two dimensional Laplace equation for the flow induced by 
the vortex sheet. 

In a finite domain R with boundary S equations (l)-(5) are equiva- 
lent to the Bateman variational principle that 

(6) I := I p dR • 

R 

is stationary- In fact according to equations (3) and (4) , a varia- 
tion causes a variation 

6p = - pq • 6q 

Thus 

iSl = - I pg*V 5«() dR 
R 

=! I (S(|) 7*(pg) dR - I 6t|) p q^ ds 
R S 

and the boundary terms vanish if 6(|) = 0 or g^ - 0. 


NUMERICAL SCHEME 

The Bateman variational principle v/ill be used to derive differ- 
ence formulas through the introduction of a discrete approximation to 
the integral I defined by equation (6) . This leads to a central diff- 
erence scheme. When such a scheme is used to compute the flow past a 
profile with fore and aft symmetry, such as an ellipse, the fore and 
aft symmetry is preserved in the solution, and expansion shocks will 
appear in transonic flow. Thus any scheme which is not desymmetrized 
in some way is restricted to subsonic flow. The basic difference 
formulas will therefore be modified by the addition of artificial vis- 
cosity to introduce the desired directional bias in the supersonic 
zone. 

In order to represent the Bateman integral, the region in which 
the flow is to be computed is divided into distorted cubic cells. 


generated from cubes by separate transformations between local 
coordinates and Cartesian coordinates x,y,z, as illustrated 

in Figure 1. 




Figure 1 


The vertices of the cells define the computational mesh, and subscripts 
i/jfk will be used to denote the value of a quantity at a mesh point. 

In order to reduce the amount of computation a simple one point inte- 
gration scheme v;ill be used, in which the contribution of each cell 
to the integral will be evaluated as the pressure at the cell center 
(defined as the point mapped from the center of the cube in the X,Y,Z 
coordinate .system) multiplied by the cell volume. Quantitites eval- 
uated at the cell centers will be denoted by subscripts i+1/2, j+1/2, 
k+1/2 . Averaging and difference operators will be introduced through 
the notation 


^X^i,j,k ' 2 ^^i+l/2,j,k"^ ^i-l/2,j,k^ 

'^X^i,j,k " ^i+l/2,j,k " ^i-l/2,j,k 

It will also be convenient to use notations such as 


xx^ *'x**'x*' ' 

^XY^ = yjr(PYf) 

xx^ = '*X<V> ' 



Numbering the vertices of a particular cell from 1 to 8 as in 


If X . ,y . , z . 

1 1 ' X 


are the Cartesian coordi- 


Figure 1, the vertices in the local coordinates are assumed to be at 

Xi - ± I , Vi = t I , Z = 1 i . _ _ _ 

nates of the i vertex, the local mapping is then defined by the 
trilinear form 


(7) 


X 


-8 I -I- x.x)(ii. ViV)li.+ z.z) 


with similar formulas for y, z. The potential is assumed to have a 
similar form inside, the cell 


( 8 ) 


4. -’s + x.xj(|i- XiV)(i+ z.x) 


These formulas preserve the continuity of x,y,i; and (}) at thp cell 
boundaries because the mappings in each cell reduce -to the same bilin- 
ear form at the conmion face. At the cell center the derivatives of 
the transformation can be evaluated by formulas such as 

i 1 . 

! ” 7 1^2 ■ =1 + ==4 ■ ==3 ^0 ■ ^ 7 ) “ "yz^x^ 

similarly it follows from equation (8) that 


' '^’XY " ’‘z^XX*^ ' ‘^'XYZ ^XYZ'^ 

In order to evaluate the contribution of each cell to the Bateman 
integral it is now necessary to express the pressure and cell volume in 
terms of the local derivatives of the mapping and the potential. Let II 
be the transformation matrix 

(9) II = 

♦ 

and let h be the determinant of H. Then the metric tensor is defined 
by the matrix 

(10) G = 

Also the contra variant velocity components are U,V,V'J where 



( 11 ) 


“ U ■ 

■ 

-1 

'♦x' 

V 

= G 

<l>y 



- ‘I'z • 


-f W<|>2 


Then 


and the variation in p duo to a variation 6»|i is 




According to the one point integration scheiuo, the contribution from 
the coll centered at i+1/2, jH-1/2, k+1/?. is the volume of the cell, 
given by the determinant j-i*l/2 k+1/2 ' multiplied by the pres- 

sure Pi^.i/2,jH.i/2,ki«l/2* 


On setting 


3X 


3 (^ 


i.j/k 


0 


and collecting the contributions from the 8 cells with a conunon vertex 
i,j,k, wo then obtain the formula 


(12) lygfijjCphU) + (phV) + Uj^j^, 62 (phW) = 0 

at each interior mesh point. Along the boundary there are only 4 
cells adjacent to each mesh point, and equation (12) is correspondingly 
modified. 

Equation (12) is a discrete approximation to i!ie conservation law 

(13) . I 37 (phU) -I* >1^ (phV) -t (phW) 

i 

which can bo derived directly from equation (1) by using the tensor 
formula for the divergence operatox* [9J . In fact wo can derive equa- 
tion (12) by representing a flux balance through a set of aujriliary 
cells, each of which is generated from a cube joining the ceiiters of 8 
primary colls, as illustrated in Pigiu-e 2. 



Figure 2 





1 


In this interpretation ^ PhtJ) ^^^2 j+1/2 )'‘}-l/2 approximation to 

the flux across the face X - 0 of that part of the secondary cell v/hich 
lies in the primary cell A in Figure 2. The boundary condition (5) 
reduces to U = 0/ V = 0 or W = 0 on cell faces which coincide with the 
boundary. The flux balance is then represented with secondary cells 
bounded on one or more faces by the body surface^ as illustrated in 
Figure 3. 


Flux balance cell 



, Figure 3 

The lumping error introduced by the one point integration scheme 
now appears as an error introduced by calculating p, h, U, V, at the 
corners of each secondary cell, instead of averaging these quantities 
over the cell faces. If the vertices of the primary cells are gener- 
ated by a global mapping smooth enough to allow Taylor series expan- 
sions of x,y,2 as functions of X,Y,Z, then the contributions to this 
error from adjacent primary cells offset each other, with the result 
that equation (12) approximates equation (13) with a second order local 
discretization error. 

The one point integration scheme has another disadvantage/ however, 
which can be seen from the following simple example. Setting h = 1, 
p — 1, equation (12) reduces in the two dimensional case to 


^'^Yy'^XX ^XX'^YY^'^ “ ° 


which is the rotated Laplacian scheme 


, '*’i+l,j+l ‘*’i-l,j+l ‘^i+i;j-i ‘^i-l,j-l “ 


= 0 



The odd nnd even points are decoupled, 30 that high Hrequoncy oscil- 
lations in v/hich «}* ” i odd points and -1 at even points are admitted 
by the scheme. To overcome this difficulty we can shift the point of 
evaluation of the flux across the side AB in Figure A from A towards 
the center of the side by adding u compensating term - 



Flux at C is 





» 


Figure 4 

The addition of similar terms on all faces produces the formula 

'"XX^YY ■ ^'^XYXY^'*' “ ° 

which reduces to the usual 5 point second order accurate formula when 
e = 1/2, and to the 9 point fourth order accurate formula when e = 1/3. 

Similarly in the discrete approximation to equation (13) we want 
to prevent excessive spatial averaging in the approximation of 
*^XX'*^YY'‘^ZZ" allowing for the dependence of p on <{>x^ coef- 
ficients of equation (13) are 

= ph(g^^ - uVa^) 

A^ = ph(g^2 - vVa^) 

Ag = ph(g^^ - 

where g^^ are the elements of G~^. In order to compensate equation 

(12), we can use these coefficients to determine the magnitude of the 

♦ 

terms which should be added to shift the locations at which <{>„ / ‘I’v f ‘J’r, 
are effectively evaluated in calculating the fluxes aci'oss each face, 
for example, e shift in the Y direction. Collecting 

the contributions from each of the 8 primary cells surrounding a 
mesh point, we obtain the following formulas. Let 



Qjjy “ (Ajj + 


with similai: formulas for , and let 


'^XYZ ° * ^'^XYZ* 


Then the final compensated equation is 


(14) g^gCj^CphU) + + yjjy 62 (phW) 


e{p2<5y^yQxY+ ^x^YZ^YS'** '^y'^ZX^ZX'’ 2 ^XYZ^XYZ^ ~ ° 


where 0 e £ 1-/2.. in practice the value e = 1/2 has been used. 

It remains to add the artificial viscosity required to desymmetrize 
the scheme iJi the supersonic zone. Instead of equation (13) we shall 
satisfy the modified conservation law 


(phU-f-P) + Iy (phV+Q) + Iy (phW+R) = 0 


where the added fluxes P/Q,R are proportional to the cell widths in 
the phj'sical domain, with the result that the correct conservation law 
is recovered in the limit as the cell width is reduced to zero. The 
artificial viscosity is designed to produce an effective switch to 
upwind differencing in the supersonic zone. Presuming the distribution 
of mesh, points to be smooth, it is constructed in the following manner. 
First ve introduce the sv/itching function 


U = h max "* I 

A 

and v?e constiruct P by the formula 


P = VI 




and Q, R by similar forraulas. Then 


^i+l/2,j,k 


i / j / 


^i+l,j,k 


if U > 0 
if U < 0 


with similar shifts for Q, R. Finally, equation (14) is modified by 
the addition of 








'i 


since M = 0 when, q < a, P,Q,R vanish in the subsonic zone. In the 
supersonic zone th4.'*y approxiraato -ii(u|<SyP / “IJlvjfiyp , -p] V?| 6,,p . 

It may bo verified [11] that the coefficients of the third deriva- i 

tives of tlie potential such as ‘•'ZXX ir^troduced by 

P/Q/R ore the same as in the artificial viscosity generated by the 

rotated difference scheme which has been previously used in three 
dimensional transonic flow calculations [1,2,4]. 

Finally the nonlinear equations generated by this discretisation 
process are solved by a generalized relaxation method v/hich is derived 
by embedding the steady state equation in an artificial time dependent 
equation. Thus we solve a discrete approximation to 

1^ (phU+P) + Iy (Phv+Q) + Iy (phV7+R) 

where the coefficients a,^,y are chosen to make the flow direction 
timelike, as in the steady state equation, and 6 controls the damp- 
ing [41. 

4. CONSTRUCTION OP THE MESH 
% 

She formulation of the artificial viscosity presupposes a smooth 
distribution of cells. Also the one point integration scheme will 
cause a loss of accuracy if the mesh is not smooth. It is important, 
therefore, to use a reasonably smooth mesh. This is most easily 
accomplished by using global mappings to generate the mesh points. 

All other steps, such as the transformation of the equations of motion, 
are then taken over by the numerical scheme. 

Swept v/ing calculations have been performed on a mesh generated 
by a sheared parabolic coordinate system, which has ber^n found to give 
good results vrith earlier methods [1,2,4]. First, we introduce 
parabolic coordinates in planes containing the wing section by the 
transformation 

_ _ • ' V2 

X -l- iY = { [x - Xq (z) + i(y - y^ (z) ) ] /t (z) } 

Z = ,z . ’ 

where z is the spanwiso coordinate, Xq{z) and yQ(2) define a singular 
line just inside the loading edge, and t(z) is a scaling factor which 


can be used to control the number of colls covering the wing. 



x,y,2 


X,Y,Z 


Figure 5 

The effect of this transformation is to unwrap tho wing to form a 
shallow bump Y = S(K,Z) , as illustrated in Figure 5. Then we use a 
sheai'ing ' ?v,.;,sformation 

X = X , y = 5 - s(x,z) , z = z 

to map the wing to the surface Y = 0. The mesh is now constructed by 
the reverse sequence of transformations from a rectangular grid in tho 
X/Y,z coordinate system. The vortex sheet trailing behind tho wing is 
assumed to coincide witli the coordinate surface leaving tlie trailing 
edge. This mesh can bo modified to treat wing cylinder combinations 
by first mapping the cylinder to a vortical slit by a Joukowsky trans- 
formation, as illustrated in Figure 6, and then using tho same sequence 
of transforp’.at.ions to generate a sheared parabolic coordinate around 
the wing projecting from the slit. 



Front view of 
wing-cylinder combination 



Cylinder mapped to vertical slit 


Figure 6 








An alternative mesh generating scheme for v;ing fuselage combina- 
tions [12] starts by introducing cylindrical coordinates as illustrated 
in Figure 7. 

0 - u/2 



Cylindrical coordinate system for 
wing body combination 

Figure 7 

In each cylindrical surface the v/ing section then appears as a profile 
in a channel bounded by the intersection of the cylinder with the 
plane of symmetry at 0 = + tt/2. This configuration can be mapped to 
a channel with a bump on the upper wall, as illustrated in Figure 8, 
by the transformation 


a ~ log [l - cosh(C)} 


0 = ir/2 




e = -ir/2 


0 = 7T/2 


0 = -tt/2 


0 


C . ' 


Figure 8 

Finally the bump is removed by a shearing transformation. 


5. 


UKSULTS 


ivo exauiplcjs o£ nuuorical calculaUionn pror.cnUod in this sec- 
tion Uo iXiustraUe the capability oil the iinite voluiae method. The 
first nosh genejrating procedure proposed in Section 4, the sheared 
parabolic coordinate system, was used in these calculations, both of 
which were performed on a sequence, of three, progressively finer meshes. 
After the. calculation on each of the first two meshes, the number of 
intervals was doubled in each coordinate direction and the interpolated 
result was used as the starting point for the calculation on the now 
mesh. The fine mesh contained 160 intervals in the chox'dwiso X direc- 
tion, 16 intervals in the normal V direction, and 32 intervals in the 
spanwise Z direction, for a total of 81,920 cells. 100 relaxation 
cycles wore used on each mesh. Such a calculation takes about 15 min- 
utes on a CDC 7600. 

The first example is a calculation of the flow past the ONRRA ll6 
wing, for v/hich experimental data is available [13], The result is 
displayed in Figure 9. Separate pressure distributions are shown for 
stations at 20, <15, 65 and 95 percent of the se.mispnn. Section lift 
and drag coefficients CL and CD wore obtained by integrating the pres- 
sure coefficient CP over the profile. The critical pressure coeffici- 
ent at which the flow has sonic speed is marked by a horizontal line 
on the pressure axis. Although the calculation did not include a 
boundary layer correction, it can be seen that the agreement with the 
oxpcriwcntal data is quite good. The triangular shock pattern is 
clearly visible in the three dimensional plot of the pressure distri- 
bution (Figure 9f ) . The front shock, emanating from the leading edge 
at the wing root, merges with the rear shock about three quarters of 
the way out across the span. The second example is indicative of the 
level of geometric complexity which can be treated with the existing 
code. The result is displayed in Figure 10. It is for a Douglas DC 10 
wing mounted on a cylinder in a low mid position. The true DC 10 coii- 
figuriition is not cxact.ly modeled, because the code does not provide 
for a viing root fillet. 

These resu'.ts confirm the promise of the now method. It appears 
that it can be used to treat configurations of more or less ai'bitrary 
complexity, subject to limits sot by the power of the available comput- 
ers. The extension to new configurations is primarily a matter of 
devising mesh generating schemes, since the internal computations are 
essentially independent of the configuration, apart from the identi- 
fication of which elements are the boundary elements. 
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